06. 행렬 함수¶
6.1 기본 연산¶
t()¶
t(x)는 행렬 x의 전치행렬(행과 열이 서로 바뀜)을 반환합니다.
인자
x: 전치할 행렬(또는 데이터프레임). 벡터를 넣으면 1행 n열의 행렬로 취급되어 전치 결과는 n행 1열이 됩니다.
(x <- matrix(1:12, 3, 4))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 4 7 10
#> [2,] 2 5 8 11
#> [3,] 3 6 9 12
t(x)
#> [,1] [,2] [,3]
#> [1,] 1 2 3
#> [2,] 4 5 6
#> [3,] 7 8 9
#> [4,] 10 11 12
diag()¶
diag(x = 1, nrow, ncol, names = TRUE)는 대각행렬을 반환합니다.
인자
x: 무엇을 넣는지에 따라 함수의 동작이 완전히 달라지는 다형(polymorphic) 인자입니다.- 숫자 하나(예:
3)를 넣으면 그 크기의 단위행렬(주대각선이 모두 1인 정방행렬)을 만듭니다. - 벡터(예:
c(2, 5, 3))를 넣으면 그 값들을 주대각선 원소로 하는 대각행렬을 만듭니다. - 행렬을 넣으면 반대로, 그 행렬의 주대각선 원소들만 뽑아 벡터로 돌려줍니다.
- 기본값
1은 인자를 아무것도 지정하지 않은diag()가 1×1 단위행렬(즉matrix(1))이 되게 합니다.
- 숫자 하나(예:
nrow,ncol: 결과 행렬의 행 수·열 수를 직접 지정합니다. 생략하면 x의 길이에 맞춰 정방행렬이 만들어집니다.names: x에 이름(names 속성)이 있을 때 그 이름을 결과 행렬의 dimnames로 유지할지 여부입니다. 기본값 TRUE.
# 주대각선의 원소가 모두 1인(단위행렬) 3x3 대각행렬
diag(3)
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
# 주대각선의 원소가 모두 10인 3x4 대각행렬
diag(10, 3, 4)
#> [,1] [,2] [,3] [,4]
#> [1,] 10 0 0 0
#> [2,] 0 10 0 0
#> [3,] 0 0 10 0
# 주대각선의 원소를 지정한 대각행렬
diag(c(2, 5, 3, 1))
#> [,1] [,2] [,3] [,4]
#> [1,] 2 0 0 0
#> [2,] 0 5 0 0
#> [3,] 0 0 3 0
#> [4,] 0 0 0 1
# 행렬 x의 대각원소만 뽑아내기
(x = matrix(c(1,3,-2, 5, 7, -3, 1, 0, 1), 3, 3))
diag(x)
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
#> [1] 1 7 1
diag(x) <- value처럼 대입 형태로 쓰면 기존 행렬의 대각원소만 골라서 일괄 수정할 수도 있습니다.
y <- matrix(0, 3, 3)
diag(y) <- c(1, 2, 3)
y
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 2 0
#> [3,] 0 0 3
upper.tri(), lower.tri()¶
upper.tri(x, diag = FALSE)와 lower.tri(x, diag = FALSE)는 행렬 x에서 대각선을 기준으로 위쪽(상삼각) 또는 아래쪽(하삼각)에 해당하는 위치를 TRUE로 표시한, x와 같은 크기의 논리 행렬을 반환합니다.
인자
x: 크기를 기준으로 삼을 행렬. 실제 값 자체는 결과에 쓰이지 않고, 오직 행 수·열 수만 참조됩니다.diag: 대각선 자신을 결과에 포함할지 여부입니다. 기본값 FALSE는 대각선을 제외하며, TRUE로 지정하면 대각선도 TRUE로 표시됩니다.
행렬 자체를 바꾸는 함수가 아니라, 조건 추출용 논리 행렬을 만들어 준다는 점에서 m3[m3 > 5] 같은 조건 인덱싱과 짝을 이룹니다. 상관행렬처럼 대각선을 기준으로 위아래가 대칭인 행렬에서 중복 없이 한쪽만 뽑아낼 때 특히 유용합니다.
(m <- matrix(1:9, 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 4 7
#> [2,] 2 5 8
#> [3,] 3 6 9
upper.tri(m)
#> [,1] [,2] [,3]
#> [1,] FALSE TRUE TRUE
#> [2,] FALSE FALSE TRUE
#> [3,] FALSE FALSE FALSE
m[upper.tri(m)]
#> [1] 4 7 8
lower.tri(m, diag = TRUE)
#> [,1] [,2] [,3]
#> [1,] TRUE FALSE FALSE
#> [2,] TRUE TRUE FALSE
#> [3,] TRUE TRUE TRUE
%*%¶
x %*% y는 행렬 x와 행렬 y의 곱셈(행렬곱) 결과를 반환합니다.
인자
x,y: 곱셈을 할 두 행렬(또는 벡터). 원소끼리 곱하는*와 달리, x의 열 수와 y의 행 수가 반드시 같아야 계산할 수 있습니다. x가 m×n 행렬이고 y가 n×p 행렬이면, 결과는 m×p 행렬이 됩니다.
(x <- 1:4)
#> [1] 1 2 3 4
(y <- diag(x))
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 0 0
#> [2,] 0 2 0 0
#> [3,] 0 0 3 0
#> [4,] 0 0 0 4
(z <- matrix(1:12, ncol = 3, nrow = 4))
#> [,1] [,2] [,3]
#> [1,] 1 5 9
#> [2,] 2 6 10
#> [3,] 3 7 11
#> [4,] 4 8 12
x %*% x
#> [,1]
#> [1,] 30
y %*% z
#> [,1] [,2] [,3]
#> [1,] 1 5 9
#> [2,] 4 12 20
#> [3,] 9 21 33
#> [4,] 16 32 48
y %*% x
#> [,1]
#> [1,] 1
#> [2,] 4
#> [3,] 9
#> [4,] 16
x %*% z
#> [,1] [,2] [,3]
#> [1,] 30 70 110
crossprod(), tcrossprod()¶
crossprod(x, y = NULL)는 t(x) %*% y를, tcrossprod(x, y = NULL)는 x %*% t(y)를 계산합니다.
인자
x: 기준이 되는 행렬y: 함께 곱할 행렬. 기본값 NULL로 두면 y에 x가 그대로 쓰여,crossprod(x)는t(x) %*% x,tcrossprod(x)는x %*% t(x)와 같아집니다.
회귀분석의 정규방정식(normal equation)처럼 \(X^{T}X\) 형태의 계산은 통계 모형 적합에서 매우 자주 등장합니다. t(x) %*% x로 직접 계산해도 결과는 같지만, crossprod(x)는 전치행렬을 별도 객체로 만들지 않고 내부적으로 한 번에 계산하므로 특히 행렬 크기가 클 때 더 빠르고 메모리도 적게 씁니다.
(a <- matrix(1:6, nrow = 2))
#> [,1] [,2] [,3]
#> [1,] 1 3 5
#> [2,] 2 4 6
crossprod(a) # t(a) %*% a 와 동일
#> [,1] [,2] [,3]
#> [1,] 5 11 17
#> [2,] 11 25 39
#> [3,] 17 39 61
t(a) %*% a
#> [,1] [,2] [,3]
#> [1,] 5 11 17
#> [2,] 11 25 39
#> [3,] 17 39 61
tcrossprod(a) # a %*% t(a) 와 동일
#> [,1] [,2]
#> [1,] 35 44
#> [2,] 44 56
a %*% t(a)
#> [,1] [,2]
#> [1,] 35 44
#> [2,] 44 56
outer()¶
outer(X, Y, FUN = "*", ...)는 두 벡터 또는 배열의 외적(outer product)을 반환합니다. 외적이란 X의 모든 원소와 Y의 모든 원소를 짝지어 FUN을 적용한 결과입니다. 기본값(곱셈)에서는 X에 Y의 전치행렬을 곱한 것과 같으며, outer(x, y)는 x %*% t(y), x %o% y와 동일합니다.
인자
X,Y: 외적을 계산할 두 벡터(또는 배열). 결과는 length(X)행 × length(Y)열의 행렬이 됩니다.FUN: X의 각 원소와 Y의 각 원소 쌍에 적용할 함수입니다. 기본값"*"(곱셈) 대신"+"같은 다른 연산자나 직접 정의한 함수(16장 참고)를 넣을 수도 있습니다....: FUN에 추가로 전달할 인자입니다.
(x <- 1:4)
#> [1] 1 2 3 4
(y <- 5:8)
#> [1] 5 6 7 8
outer(x, y)
#> [,1] [,2] [,3] [,4]
#> [1,] 5 6 7 8
#> [2,] 10 12 14 16
#> [3,] 15 18 21 24
#> [4,] 20 24 28 32
x %*% t(y)
#> [,1] [,2] [,3] [,4]
#> [1,] 5 6 7 8
#> [2,] 10 12 14 16
#> [3,] 15 18 21 24
#> [4,] 20 24 28 32
x%o%y
#> [,1] [,2] [,3] [,4]
#> [1,] 5 6 7 8
#> [2,] 10 12 14 16
#> [3,] 15 18 21 24
#> [4,] 20 24 28 32
FUN을 기본값인 *에서 +로 바꾸면 곱셈 대신 덧셈으로 조합됩니다.
outer(x, y, FUN = "+")
#> [,1] [,2] [,3] [,4]
#> [1,] 6 7 8 9
#> [2,] 7 8 9 10
#> [3,] 8 9 10 11
#> [4,] 9 10 11 12
kronecker() (%x%)¶
kronecker(X, Y, FUN = "*", make.dimnames = FALSE, ...)는 크로네커 곱(Kronecker product)을 반환하며, X %x% Y로도 쓸 수 있습니다.
인자
X,Y: 크로네커 곱을 계산할 두 행렬(또는 배열)FUN: 각 블록을 만들 때 사용할 연산입니다. 기본값은 곱셈("*")make.dimnames: 결과 행렬에 X와 Y의 이름을 조합한 dimnames를 붙일지 여부. 기본값 FALSE...: FUN에 추가로 전달할 인자
outer()가 두 벡터의 모든 조합을 곱해 하나의 행렬을 만드는 것이라면, kronecker()는 행렬 X의 원소 하나하나를 행렬 Y 전체로 치환하여 훨씬 큰 블록 행렬을 만듭니다. 예를 들어 X가 2×2 행렬이고 Y도 2×2 행렬이면, 결과는 4×4 행렬이 됩니다. 통계학에서는 반복측정 자료의 분산·공분산 구조를 표현하거나, 여러 그룹에 동일한 설계행렬을 반복 적용할 때 등장합니다.
(A <- matrix(1:4, 2, 2))
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 2 4
(B <- diag(2))
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
A %x% B
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 3 0
#> [2,] 0 1 0 3
#> [3,] 2 0 4 0
#> [4,] 0 2 0 4
kronecker(A, B)
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 3 0
#> [2,] 0 1 0 3
#> [3,] 2 0 4 0
#> [4,] 0 2 0 4
solve()¶
solve(a, b, tol = .Machine$double.eps, LINPACK = FALSE, ...)는 수식 a %*% x = b에서 x를 구하여 반환합니다.
인자
a: 계수행렬(정방행렬)b: 방정식 우변의 벡터 또는 행렬입니다. 생략하면 단위행렬을 넣은 것으로 간주되어, 결과적으로 a의 역행렬을 구하는 것과 같아집니다.tol: a가 특이행렬(역행렬이 존재하지 않는 행렬)에 얼마나 가까운지를 판정하는 수치 허용오차입니다.LINPACK: 더 이상 쓰이지 않는(defunct) 인자입니다. 예전 R 버전은 LINPACK이라는 계산 라이브러리를 선택할 수 있었지만, 지금은 항상 LAPACK을 사용하며 이 인자를 TRUE로 지정하면 오히려 오류가 발생합니다. 신경 쓰지 않고 생략하는 것이 좋습니다....: 그 밖에 메서드별로 추가되는 인자
(a <- matrix(c(1, 3, -2, 5, 7, -3, 1, 0, 1), 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
(b <- c(2, -1, 1))
#> [1] 2 -1 1
# ax=b에서 x를 구하여 반환
solve(a, b)
#> [1] -5 2 -3
# a의 역행렬을 반환
solve(a)
#> [,1] [,2] [,3]
#> [1,] -2.333333 2.666667 2.333333
#> [2,] 1.000000 -1.000000 -1.000000
#> [3,] -1.666667 2.333333 2.666667
참고:
solve()는 내부적으로 역행렬을 직접 구하지 않고 LU 분해를 이용해 연립방정식을 풀며, a가 정방행렬이 아니거나 특이행렬이면 오류가 발생합니다.
backsolve(), forwardsolve()¶
backsolve(r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)와 forwardsolve(l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)는 각각 상삼각행렬 r, 하삼각행렬 l에 대해 r %*% z = x(또는 l %*% z = x)를 만족하는 z를 구합니다.
인자
r/l: 삼각행렬.backsolve()는 기본적으로 상삼각행렬을,forwardsolve()는 하삼각행렬을 기대합니다.x: 방정식 우변의 벡터 또는 행렬k: r(또는 l)에서 실제로 사용할 행·열의 개수입니다. 기본값은 전체 열 수이며, 일부만 쓰고 싶을 때 줄일 수 있습니다.upper.tri: 입력 행렬을 상삼각으로 볼지(TRUE) 하삼각으로 볼지(FALSE) 지정합니다.backsolve()는 기본값이 TRUE,forwardsolve()는 기본값이 FALSE로 서로 반대입니다.transpose: TRUE로 지정하면 r(또는 l)을 전치한 것으로 간주하고 계산합니다.
solve()로도 같은 문제를 풀 수 있지만, 삼각행렬은 위(또는 아래)에서부터 한 줄씩 대입만 하면 답이 바로 나오는 특수한 구조를 가지고 있습니다. backsolve()·forwardsolve()는 이 구조를 활용해 일반적인 solve()보다 훨씬 적은 연산으로 더 빠르게 계산합니다. qr()이나 chol()로 분해한 결과를 뒤이어 처리할 때 자주 함께 쓰입니다.
# 상삼각행렬
(U <- matrix(c(2, 0, 1, 3), nrow = 2))
#> [,1] [,2]
#> [1,] 2 1
#> [2,] 0 3
backsolve(U, c(14, 18))
#> [1] 4 6
# 하삼각행렬(U의 전치)
(L <- t(U))
#> [,1] [,2]
#> [1,] 2 0
#> [2,] 1 3
forwardsolve(L, c(8, 22))
#> [1] 4 6
6.2 행·열 통계¶
행렬이나 데이터프레임에서 행별·열별 합계나 평균을 구할 때, for() 반복문을 직접 돌리거나 apply(x, 1, sum)처럼 apply 계열 함수를 쓸 수도 있습니다. 하지만 아래 함수들은 이 계산 전용으로 내부가 최적화되어 있어 더 빠르고 코드도 간결합니다.
rowSums()¶
rowSums(x, na.rm = FALSE, dims = 1)는 숫자를 담고 있는 2차원 이상의 배열 또는 데이터프레임 x의 행 합계를 반환합니다.
인자
x: 행렬, 배열, 또는 데이터프레임na.rm: 결측값(NA)을 계산에서 제외할지 여부입니다. 기본값 FALSE에서는 해당 행에 NA가 하나라도 있으면 합계 결과도 NA가 됩니다.dims: 몇 번째 차원까지를 묶어 "행"으로 취급할지 지정합니다. 일반적인 2차원 행렬에서는 기본값 1로 충분하며, 3차원 이상 배열을 다룰 때 의미가 커집니다.
(x <- matrix(c(1, 3, -2, 5, 7, -3, 1, 0, 1), 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
rowSums(x)
#> [1] 7 10 -4
a <- 1:5
b <- 1:5*2
c <- 1:5*3
df <- data.frame(a, b, c)
df
#> a b c
#> 1 1 2 3
#> 2 2 4 6
#> 3 3 6 9
#> 4 4 8 12
#> 5 5 10 15
rowSums(df)
#> [1] 6 12 18 24 30
apply(x, 1, sum) 방식은 rowSums() 함수보다 더 느리지만 sum 대신 다른 함수도 자유롭게 쓸 수 있습니다.
colSums()¶
colSums(x, na.rm = FALSE, dims = 1)은 숫자를 담고 있는 2차원 이상의 배열 또는 데이터프레임 x의 열 합계를 반환합니다. 인자의 의미는 rowSums()와 동일하며, 방향만 열 기준으로 바뀝니다.
(x <- matrix(c(1, 3, -2, 5, 7, -3, 1, 0, 1), 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
colSums(x)
#> [1] 2 9 2
a <- 1:5
b <- 1:5*2
c <- 1:5*3
df <- data.frame(a, b, c)
df
#> a b c
#> 1 1 2 3
#> 2 2 4 6
#> 3 3 6 9
#> 4 4 8 12
#> 5 5 10 15
colSums(df)
#> a b c
#> 15 30 45
주의 —
rowsum()과 헷갈리지 않기: 대소문자만 다른rowsum(x, group, reorder = TRUE, ...)이라는 함수가 별도로 존재합니다.rowSums()는 각 행의 합계를 그대로 구하는 반면,rowsum()은group인자로 지정한 그룹(범주)별로 행들을 묶어 합산합니다.group은 각 행이 어느 그룹에 속하는지 나타내는 벡터이며,reorder(기본값 TRUE)는 결과 행을 그룹 이름 순서대로 재배열할지를 정합니다. 이름이 비슷해 혼동하기 쉬우므로 원하는 계산이 "행 전체의 합"인지 "그룹별 합"인지에 따라 정확히 구분해서 써야 합니다.
rowMeans()¶
rowMeans(x, na.rm = FALSE, dims = 1)는 숫자를 담고 있는 2차원 이상의 배열 또는 데이터프레임 x의 행 평균을 반환합니다. 인자의 의미는 rowSums()와 같습니다.
(x <- matrix(c(1, 3, -2, 5, 7, -3, 1, 0, 1), 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
rowMeans(x)
#> [1] 2.333333 3.333333 -1.333333
a <- 1:5
b <- 1:5*2
c <- 1:5*3
df <- data.frame(a, b, c)
df
#> a b c
#> 1 1 2 3
#> 2 2 4 6
#> 3 3 6 9
#> 4 4 8 12
#> 5 5 10 15
rowMeans(df)
#> [1] 2 4 6 8 10
colMeans()¶
colMeans(x, na.rm = FALSE, dims = 1)은 숫자를 담고 있는 2차원 이상의 배열 또는 데이터프레임 x의 열 평균을 반환합니다. 인자의 의미는 rowSums()와 같으며, 방향만 열 기준입니다.
(x <- matrix(c(1, 3, -2, 5, 7, -3, 1, 0, 1), 3, 3))
#> [,1] [,2] [,3]
#> [1,] 1 5 1
#> [2,] 3 7 0
#> [3,] -2 -3 1
colMeans(x)
#> [1] 0.6666667 3.0000000 0.6666667
a <- 1:5
b <- 1:5*2
c <- 1:5*3
df <- data.frame(a, b, c)
df
#> a b c
#> 1 1 2 3
#> 2 2 4 6
#> 3 3 6 9
#> 4 4 8 12
#> 5 5 10 15
colMeans(df)
#> a b c
#> 3 6 9
6.3 차원·분해¶
nrow(), ncol()¶
nrow(x)는 행의 개수를 반환합니다.ncol(x)는 열의 개수를 반환합니다.NROW(x)는 행의 개수를 반환합니다. (벡터 계산 가능)NCOL(x)는 열의 개수를 반환합니다. (벡터 계산 가능)
인자
x: 행 수·열 수를 알고 싶은 행렬, 배열, 데이터프레임 등. 인자는 이것 하나뿐이지만, x의 종류(행렬인지 벡터인지, 몇 차원 배열인지)에 따라 결과가 달라지므로 아래에서 사례별로 살펴봅니다.
(x <- array(1:24, dim = c(3,4,2))) # 3행*4열의 3차원 배열
#> , , 1
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 1 4 7 10
#> [2,] 2 5 8 11
#> [3,] 3 6 9 12
#>
#> , , 2
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 13 16 19 22
#> [2,] 14 17 20 23
#> [3,] 15 18 21 24
#>
nrow(x)
#> [1] 3
ncol(x)
#> [1] 4
3차원 이상 배열에서도 nrow()·ncol()은 항상 첫 번째·두 번째 차원(행·열)만 알려 줄 뿐, 세 번째 차원(위 예제에서는 2) 이상은 별도로 dim(x)[3]처럼 확인해야 합니다.
NROW()·NCOL()은 행렬이 아닌 벡터도 계산할 수 있다는 점이 nrow()·ncol()과 다릅니다. 벡터는 "행"이라는 개념이 없으므로 nrow()는 NULL을 반환하지만, NROW()는 벡터를 원소가 하나뿐인 열로 취급해 길이를 그대로 돌려줍니다.
det()¶
det(x, ...)는 행렬 x의 행렬식을 구하여 반환합니다. 행렬식은 정방행렬(n×n 행렬)인 경우에만 구할 수 있습니다.
인자
x: 행렬식을 구할 정방행렬...: 예전 버전에는 계산 방식을 고르는method인자가 있었지만 현재는 쓰이지 않으며, 값을 넣어도 오류 없이 무시됩니다.
행렬식이 0이면 그 행렬은 역행렬이 존재하지 않는 특이행렬(singular matrix)이며, solve()로 역행렬을 구하려 하면 오류가 발생합니다.
determinant()¶
determinant(x, logarithm = TRUE, ...)는 행렬식을 부호와 크기로 나누어 반환합니다. 내부적으로는 LU 분해를 이용해 계산하며, det()는 사실 이 함수를 호출한 뒤 부호와 크기를 다시 곱해 값 하나로 합쳐 주는 간단한 래퍼(wrapper)입니다.
인자
x: 행렬식을 구할 정방행렬logarithm: TRUE(기본값)이면 행렬식 절댓값에 로그를 취한 값을(modulus), FALSE이면 로그를 취하지 않은 절댓값 그대로를 반환합니다. 부호(양수/음수)는logarithm설정과 무관하게 항상sign에 별도로 담겨 반환됩니다....: 추가 메서드별 인자
행렬 크기가 커지면 행렬식의 실제 값이 지나치게 크거나 작아져(오버플로·언더플로) det()로는 정확한 값을 구하기 어려운 경우가 생깁니다. determinant()는 부호와 절댓값의 로그를 따로 반환하므로, 이런 극단적인 값도 안정적으로 다룰 수 있습니다.
(m <- matrix(c(4, 2, 7, 6), 2, 2))
#> [,1] [,2]
#> [1,] 4 7
#> [2,] 2 6
det(m)
#> [1] 10
determinant(m, logarithm = TRUE)
#> $modulus
#> [1] 2.302585
#> attr(,"logarithm")
#> [1] TRUE
#>
#> $sign
#> [1] 1
#>
#> attr(,"class")
#> [1] "det"
# 로그값을 다시 지수변환하면 det()와 같은 값이 나옴
exp(determinant(m, logarithm = TRUE)$modulus)
#> [1] 10
#> attr(,"logarithm")
#> TRUE
eigen()¶
eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE)는 행렬 x의 고유값과 고유벡터를 반환합니다.
고유값과 고유벡터는 행렬 A에 대하여 \(Av = \lambda v\) 등식을 만족하는 상수와 열벡터입니다. 등식에서 v는 고유벡터이고 \(\lambda\)는 고유값입니다. 고유값과 고유벡터는 정방행렬(n×n 행렬)인 경우에만 구할 수 있습니다. 데이터분석에서는 공분산행렬을 고유값 분해하여 주성분분석(PCA)의 주성분을 구하는 등 차원축소에 널리 쓰입니다.
인자
x: 고유값·고유벡터를 구할 정방행렬symmetric: x가 대칭행렬인지를 지정합니다. 생략하면 R이 x가 대칭인지 자동으로 판단하지만, x가 대칭임을 미리 알고 있다면 직접 TRUE로 지정하는 편이 더 안정적인 전용 알고리즘을 쓰게 되어 계산도 더 빠르고 정확합니다.only.values: TRUE로 지정하면 고유벡터는 계산하지 않고 고유값만 반환합니다. 고유벡터가 필요 없을 때 계산량을 줄일 수 있습니다.EISPACK: 더 이상 쓰이지 않는(defunct) 인자로, 값을 지정해도 무시됩니다. 예전에는 EISPACK이라는 계산 라이브러리를 선택하는 용도였지만, 지금은 항상 더 최신인 LAPACK을 사용합니다.
(x <- matrix(1:4, nrow = 2))
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 2 4
eigen(x)
#> eigen() decomposition
#> $values
#> [1] 5.3722813 -0.3722813
#>
#> $vectors
#> [,1] [,2]
#> [1,] -0.5657675 -0.9093767
#> [2,] -0.8245648 0.4159736
svd()¶
svd(x, nu = min(n, p), nv = min(n, p), LINPACK = FALSE)는 행렬 x의 특이값을 분해하여 그 결과를 반환합니다.
고유값 분해는 정방행렬(n×n 행렬)인 경우에만 가능합니다. 이를 직사각행렬(m×n 행렬)에서도 가능하게 만든 것이 특이값 분해(SVD: Singular Value Decomposition)입니다.
인자
x: 특이값 분해를 할 행렬입니다. n×p 행렬로, 정방행렬이 아니어도 됩니다.nu: 왼쪽 특이벡터(결과의$u, 행렬 U의 열)를 몇 개까지 계산할지 지정합니다. 기본값은min(n, p)로, 의미 있는 특이벡터를 모두 계산합니다.nv: 오른쪽 특이벡터(결과의$v, 행렬 V의 열)를 몇 개까지 계산할지 지정합니다. 기본값은 nu와 같습니다.LINPACK: 더 이상 쓰이지 않는(defunct) 인자입니다. TRUE로 지정하면 오류가 발생하므로, 지정하지 말고 생략하는 것이 좋습니다.
특이값을 크기순으로 늘어놓았을 때 앞쪽 몇 개만 남기고 나머지를 버리면, 원래 행렬을 가장 적은 정보 손실로 압축한 저차원 근사 행렬을 얻을 수 있습니다. 이 원리는 이미지 압축이나 추천 시스템, 잡음 제거 등에 폭넓게 활용됩니다.
x <- matrix(c(0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 0),
byrow = TRUE, nrow = 10)
x_svd <- svd(x)
x_svd$d
#> [1] 6.058225e+00 2.752981e+00 8.479424e-01 2.837992e-16 2.460311e-31
#> [6] 2.867106e-33 4.206539e-50 6.039616e-67
# 특이값 벡터 d의 첫번째와 두번째 값만 이용해 원래 행렬을 근사하기
round(x_svd$u[,c(1,2)] %*% diag(x_svd$d[c(1,2)]) %*% t(x_svd$v[,c(1,2)]))
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0 0 0 1 1 0 0 0
#> [2,] 0 0 0 1 1 0 0 0
#> [3,] 0 0 0 1 1 0 0 0
#> [4,] 1 1 1 1 1 1 1 1
#> [5,] 1 1 1 1 1 1 1 1
#> [6,] 1 1 1 1 1 1 1 1
#> [7,] 1 1 1 1 1 1 1 1
#> [8,] 0 0 0 1 1 0 0 0
#> [9,] 0 0 0 1 1 0 0 0
#> [10,] 0 0 0 1 1 0 0 0
원래 x의 2행 1열은 사실 1이었지만(다른 행과 달리 왼쪽 위 모서리가 살짝 어긋난 "잡음"이었습니다), 상위 2개의 특이값만으로 복원한 근사 행렬에서는 0으로 매끄럽게 처리되었습니다. 이는 SVD로 저차원 근사를 하면 전체적인 패턴에서 벗어난 작은 잡음이 자연스럽게 걸러지는 효과가 있음을 보여 줍니다.
qr()¶
qr(x, tol = 1e-07, LAPACK = FALSE, ...)는 행렬 x의 QR 분해 결과를 반환합니다.
QR 분해는 행렬 A를 \(A = QR\)로 분해하는 것입니다. Q는 단위 노름 직교 벡터를 갖는 행렬이고, R은 상삼각행렬입니다. 회귀분석에서 최소제곱해를 구할 때 역행렬을 직접 계산하는 것보다 수치적으로 더 안정적인 방법으로 널리 쓰입니다.
인자
x: QR 분해를 할 행렬. 반드시 정방행렬일 필요는 없습니다.tol: x의 열들 사이에 선형종속(공선성) 관계가 있는지 판정하는 허용오차입니다.LAPACK = FALSE일 때만 사용됩니다.LAPACK: TRUE로 지정하면 LAPACK 기반의 다른 알고리즘을 사용합니다. 기본값 FALSE는 R 고유의 알고리즘을 사용하며, 이때만tol에 의한 열별 선형종속 판정이 이루어집니다....: 추가 메서드별 인자
x <- matrix(c(1,2,3, 2,4,6, 3,3,3), nrow=3)
x_qr <- qr(x)
x_qr
#> $qr
#> [,1] [,2] [,3]
#> [1,] -3.7416574 -4.810702 -7.483315
#> [2,] 0.5345225 1.963961 0.000000
#> [3,] 0.8017837 0.988693 0.000000
#>
#> $rank
#> [1] 2
#>
#> $qraux
#> [1] 1.267261 1.149954 0.000000
#>
#> $pivot
#> [1] 1 3 2
#>
#> attr(,"class")
#> [1] "qr"
$rank가 2로 나온 것은 x의 두 번째 열이 첫 번째 열의 2배(선형종속)이기 때문입니다. 즉 x는 3×3 행렬이지만 실질적인 독립 정보는 2개 열 분량뿐이라는 뜻이며, 이런 상황을 QR 분해가 rank로 곧바로 알려 줍니다.
chol()¶
chol(x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...)는 대칭이면서 양의 정부호(positive definite)인 행렬 x를 숄레스키 분해(Cholesky decomposition)하여 상삼각행렬 R을 반환합니다.
인자
x: 숄레스키 분해를 할 대칭 양의 정부호 행렬입니다. 이 조건을 만족하지 않으면(대칭이 아니거나 고유값 중 하나라도 0 이하이면) 오류가 발생합니다.pivot: TRUE로 지정하면 피벗팅을 허용하여, 경계선에 있는 행렬(양의 준정부호, 고유값에 정확히 0이 섞인 경우)도 분해할 수 있습니다. 다만 이 경우 결과 R에 대해t(R) %*% R이 x와 순서가 어긋나므로,attr(R, "pivot")으로 원래 순서를 다시 맞춰 주어야 합니다.LINPACK: 더 이상 쓰이지 않는(defunct) 인자로, 지정해도 무시됩니다.tol:pivot = TRUE일 때 사용하는 수치 허용오차입니다....: 추가 메서드별 인자
공분산행렬은 대칭·양의 정부호 조건을 만족하는 대표적인 예입니다. eigen()이나 svd()로도 비슷한 정보를 얻을 수 있지만, 이 조건을 만족하는 행렬에 한해서는 chol()이 계산량이 훨씬 적어 다변량 정규분포 난수 생성이나 최소제곱 계산에 자주 쓰입니다.
(A <- matrix(c(4, 2, 2, 3), 2, 2))
#> [,1] [,2]
#> [1,] 4 2
#> [2,] 2 3
(R <- chol(A))
#> [,1] [,2]
#> [1,] 2 1.000000
#> [2,] 0 1.414214
# 검산: t(R) %*% R이 원래 A와 같아야 함
t(R) %*% R
#> [,1] [,2]
#> [1,] 4 2
#> [2,] 2 3
주의:
pivot = FALSE(기본값) 상태에서 x가 양의 정부호가 아니면chol()은 오류를 발생시킵니다. 공분산행렬처럼 조건이 수학적으로 보장된 경우가 아니라면, 먼저eigen(x)$values로 고유값이 모두 양수인지 확인해 보는 것이 안전합니다.
scale()¶
scale(x, center = TRUE, scale = TRUE)는 행렬이나 벡터 x를 정규화(표준화)한 결과를 반환합니다.
인자
x: 표준화할 벡터 또는 행렬center: TRUE(기본값)이면 각 열(x가 벡터면 전체)의 평균을 뺍니다. 숫자 벡터를 직접 지정하면 그 값을 빼며, FALSE로 지정하면 아예 빼지 않습니다.scale: TRUE(기본값)이면 각 열(또는 전체)의 표준편차로 나눕니다. 숫자 벡터를 직접 지정하면 그 값으로 나누며, FALSE로 지정하면 나누지 않습니다.
x가 벡터일 경우 \(z = (x-\mu)/\sigma\)를 구해서 반환합니다.
# 1~9 벡터를 표준화하기 (평균을 빼고 표준편차로 나눔)
(x <- 1:9)
#> [1] 1 2 3 4 5 6 7 8 9
scale(x)
#> [,1]
#> [1,] -1.4605935
#> [2,] -1.0954451
#> [3,] -0.7302967
#> [4,] -0.3651484
#> [5,] 0.0000000
#> [6,] 0.3651484
#> [7,] 0.7302967
#> [8,] 1.0954451
#> [9,] 1.4605935
#> attr(,"scaled:center")
#> [1] 5
#> attr(,"scaled:scale")
#> [1] 2.738613
x가 행렬일 경우에는 열을 기준으로 평균과 표준편차를 구한 후, 열 기준으로 표준화를 합니다.
(x <- matrix(1:9, ncol = 3))
#> [,1] [,2] [,3]
#> [1,] 1 4 7
#> [2,] 2 5 8
#> [3,] 3 6 9
scale(x)
#> [,1] [,2] [,3]
#> [1,] -1 -1 -1
#> [2,] 0 0 0
#> [3,] 1 1 1
#> attr(,"scaled:center")
#> [1] 2 5 8
#> attr(,"scaled:scale")
#> [1] 1 1 1
norm(), rcond()¶
norm(x, type = c("O", "I", "F", "M", "2"))는 행렬(또는 벡터) x의 크기를 나타내는 노름(norm)을 반환하고, rcond(x, norm = c("O", "I", "1"), triangular = FALSE, ...)는 행렬 x의 조건수(condition number)의 역수를 반환합니다.
인자 — norm()
x: 노름을 구할 행렬(또는 벡터)type: 계산 방식을 지정합니다."O"(또는"1", 기본값): 열 합의 절댓값 중 최댓값."I": 행 합의 절댓값 중 최댓값."F": 프로베니우스 노름(전체 원소 제곱합의 제곱근)."M": 전체 원소 중 절댓값이 가장 큰 값."2": 스펙트럴 노름(가장 큰 특이값과 같음).
인자 — rcond()
x: 조건수를 구할 정방행렬norm: 조건수 계산에 사용할 노름 종류입니다.triangular: x가 이미 삼각행렬임을 알고 있다면 TRUE로 지정해 계산을 더 빠르게 할 수 있습니다....: 추가 인자
rcond()의 값은 0에 가까울수록 x가 특이행렬(역행렬이 존재하지 않는 행렬)에 가깝다는 뜻이며, solve(x)로 역행렬을 구하는 계산이 수치적으로 불안정해질 위험이 크다는 신호입니다.
(m <- matrix(c(1, 2, 3, 4), 2, 2))
#> [,1] [,2]
#> [1,] 1 3
#> [2,] 2 4
norm(m, type = "F")
#> [1] 5.477226
norm(m, type = "1")
#> [1] 7
norm(m, type = "I")
#> [1] 6
norm(m, type = "2")
#> [1] 5.464986
# "2" 노름은 가장 큰 특이값과 정확히 같음
svd(m)$d[1]
#> [1] 5.464986
rcond() 값이 0에 극단적으로 가깝다면(예: 1e-10 이하) solve()의 결과를 그대로 믿기보다, 데이터에 중복되거나 거의 선형종속인 열이 없는지 먼저 점검해 보는 것이 좋습니다.